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Abstract 

SecDec is a program which can be used for the factorization of dimensionally 
regulated poles from parametric integrals, in particular multi-loop integrals, and the 
subsequent numerical evaluation of the finite coefficients. Here we present version 3.0 
of the program, which has major improvements compared to version 2: it is faster, 
contains new decomposition strategies, an improved user interface and various other 
new features which extend the range of applicability. 
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PROGRAM SUMMARY 


Program Title: SecDec 3.0 

Journal Reference: Comput. Phys. Comm. 196 (2015) 470. 

Catalogue identifier: AEIR_v3_0 
Licensing provisions: Standard CPC license 

Programming language: Wolfram Mathematica, perl, Fortran/C++ 

Computer: from a single PC to a cluster, depending on the problem 

Operating system: Unix, Linux 

RAM: depending on the complexity of the problem 

Keywords: Perturbation theory, Feynman diagrams, multi-loop, numerical integra¬ 
tion 

Classification: 

4.4 Feynman diagrams, 5 Computer Algebra, 11.1 General, High Energy Physics 
and Computing. 

Journal reference of previous version: Comput. Phys. Commun. 184 (2013) 2552. 
Nature of the problem: 

Extraction of ultraviolet and infrared singularities from parametric integrals ap¬ 
pearing in higher order perturbative calculations in gauge theories. Numerical inte¬ 
gration in the presence of integrable singularities (e.g. kinematic thresholds). 
Solution method: 

Algebraic extraction of singularities within dimensional regularization using iterated 
sector decomposition. This leads to a Laurent series in the dimensional regulariza¬ 
tion parameter e, where the coefficients are finite integrals over the unit-hypercube. 
Those integrals are evaluated numerically by Monte Carlo integration. The inte¬ 
grable singularities are handled by choosing a suitable integration contour in the 
complex plane, in an automated way. 

Restrictions: Depending on the complexity of the problem, limited by memory and 
CPU time. 

Running time: 

Between a few seconds and several hours, depending on the complexity of the prob¬ 
lem. 
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1 Introduction 


After the very successful Run I of the LHC, pushing the precision frontier is 
and will be one of the primary goals for the next phase of LHC data taking 
and at future colliders. 

This means that higher-order corrections within the Standard Model, and 
promising extensions beyond, need to be evaluated. In addition, in order to 
scrutinize the Higgs properties, in particular the Yukawa couplings to fermions, 
heavy quark masses should be taken into account without resorting to low 
energy approximations. Therefore, it is of primary importance to have tools 
at hand for the calculation of two- (and more) loop integrals involving several 
mass scales. 


In general, high precision calculations have in common that they involve multi¬ 
dimensional integrations over some parameters: Feynman (or Schwinger) pa¬ 
rameters in the case of (multi-)loop integrals, or parameters related to the 
integration of subtraction terms over a factorised phase space in the case 
of infrared-divergent real radiation. Usually, these calculations are performed 
within the framework of dimensional regularization, and one of the challenges 
is to factorise the poles in the regulator e. 


The program SecDec is designed to perform this task in an automated 

way, and to integrate the coefficients of the resulting Laurent series in e numer¬ 
ically, based on the sector decomposition algorithm described in Refs. [^[^, 
which was inspired by earlier ideas as contained in Refs. [^[^. Other public 
implementations of sector decomposition can be found in Refs. [8 12 . 


The numerical integration in SecDec-1.0 was restricted to Euclidean kine¬ 
matics for integrals with more than one kinematic scale. However, this restric¬ 
tion was lifted in SecDec-2.0 [2,13 , by combining sector decomposition with 


a method to deform the multi-dimensional integration contour into the com¬ 


plex plane 14 . While such a method of contour deformation already had been 
applied in various contexts - for examples at one loop see e.g. Refs. 15 21 , for 


two-loop examples Refs. 22-24 - it has been combined with the automated 


setup for the resolution of singularities and made publicly available for the 
hrst time in Ref. [^. 

Very recently, another strategy to achieve a resolution of dimensionally reg¬ 
ulated singularities in multi-loop integrals has been proposed 
utilizes integration by parts and dimension shifts. 


25,26 , which 


The original sector decomposition algorithm described in Ref. [^ is based 
on an iterative procedure, which can run into an inhnite recursion in certain 
(rather rare) situations. It was soon noticed [^ however that the structure of 
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Feynman integrals is such that a decomposition algorithm must exist which is 
guaranteed to stop, as the procedure can be mapped to a well known problem 
in convex geometry. Related observations have been made in Ref. 27 . In 
Ref. 28 , an algorithm was presented which cannot lead to inhnite recursion 


and is more efficient than previous algorithms with this property. 


In this paper, we present the implementation of two new decomposition strate¬ 
gies as alternatives to the already implemented heuristic algorithm. They are 
based on the method of Ref. and therefore guaranteed to avoid inhnite re¬ 
cursion. This is described in Section together with other new features of the 
program. As the user interface has been restructured in version 3 of SecDec, 
we describe its usage in detail in Section and give various examples in Sec¬ 
tion After the conclusions, an appendix provides a detailed description of 
the various input options and run modes. 


2 Theoretical framework 


In this section we give the expressions for the representation of multi-loop 
integrals in terms of Feynman parameters and explain our notation. 


2.1 Feynman integrals 


In order to dehne our conventions we choose a scalar integral. Integrals with 
loop momenta in the numerator only lead to an additional function of the 
Feynman parameters and invariants in the numerator and will be discussed in 
Section 13.2.41 

A scalar Feynman integral G'm. D dimensions at L loops with N propagators, 
where the propagators can have arbitrary, not necessarily integer powers Vj, 
has the following representation in momentum space: 

^- ( 1 ) 

nPpm.ipj.mj) 

J=1 

d^Ki = d^ki , Pj{{k}, {p}, m^j) = gj - i6 , 

iiT 2 

where the qj are linear combinations of external momenta pi and loop momenta 
ki. The numerical result produced by SecDec will correspond to G with 
/i = 1, except if the user specihes a different prefactor when dehning the 
graph to be calculated. 
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Denominator factors which only have a linear, but not a quadratic dependence 
on the loop momenta can appear in various contexts within loop calculations. 
In the new version presented here, SecDec can deal with such propagators. 
However, care has to be taken with the analytic continuation in this case, as 
will be discussed in Section 13.2.51 

Introducing Feynman parameters in Eq. ([^ leads to 


G = 


r(iv. 




N N 

dxj 5{l - ^Xi 




2 = 1 


■ j d^Ki... d^Ki 


n -Nu 


Y, kj M,jkj-2Yk] ■Qj + J + i5 


151 

°? N 


nf=ir(i/,) 


r(iV, - LD/2) 


i=i 


( 2 ) 


I 

0 


I'j — l 


N j/N^-{L+l)D/2 

5(1 

i=i 


jrN^-LD/2 ’ 


( 3 ) 


where 

J^(f) = det(M) 


L 


YQ,MYQ,-J-i5 

3 , 1=1 


N 

U{x) = det(M), Ny = 'Y ■ 

i=i 


( 4 ) 

( 5 ) 


In the expressions above, M is an L x L matrix containing Feynman param¬ 
eters, Q is an L-dimensional vector, where each entry is a linear combination 
of external momenta and Feynman parameters, and J is a scalar expression 
containing kinematic invariants and Feynman parameters. For more details 
and examples, see e.g. Ref. |^. 


The functions U and called graph polynomials in the following, can also be 
constructed from the topology of the corresponding Feynman graph [5,29,30 


In SecDec, this is implemented as one way to hnd the two graph polynomials. 


For a diagram with massless propagators, none of the Feynman parameters 
occurs quadratically in the function = J^q. \i massive internal lines are 

present, gets an additional term J^{x) = J^o{x) + U{x) Y, Xj-rri^,. 

i=i 

W is a positive semi-dehnite function. A vanishing U function is related to the 
UV sub divergences of the graph, a vanishing Y function to its IR divergences. 
In the region where all invariants formed from external momenta are negative, 
which we will call the Euclidean region in the following, Y is also a positive 
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semi-definite function of the Feynman parameters Xj. If some of the invariants 
are zero, for example if some of the external momenta are light-like, an IR 
divergence may appear due to vanishing. Therefore it depends on the kine¬ 
matics^ and not only on the topology (like in the UV case), whether a zero of 
T leads to a divergence or not. The necessary (but not sufficient) conditions 
for a divergence are given by the Landau equations [^[^. If all kinematic 
invariants formed by external momenta are negative, the necessary condition 
= 0 for an IR divergence can only be fulhlled if some of the parameters 
Xi are zero. These endpoint singularities can be regulated by dimensional reg- 
ularisation and factored out of the function T using sector decomposition. 
The same holds for dimensionally regulated UV singularities contained in lA. 
However, after the UV and IR singularities have been extracted as poles in 
1/e, for non-Euclidean kinematics integrable singularities related to kinematic 
thresholds remain. These singularities imply that iF is vanishing inside the 
integration region for some combinations of Feynman parameter values and 
values of the kinematic invariants. However, the integrals can be evaluated 
by deforming the integration contour into the complex plane, as explained in 
detail in Ref. [^. 


result given by SecDec will be the one for the integral as de- 


The numerica 
hned in Eq. ( ll). This implies that the prefactor 




V{Nj, — LD 12) coming 


from the Feynman parametrisation by default will be included in our numer¬ 
ical result. However, the user can dehne a different prefactor to be factored 
out from the numerical result. 


3 Structure and new features of SecDec version 3.0 


3.1 Program structure 


The program consists of two main parts, one designed for loop integrals, the 
other for more general parametric integrals. The sector decomposition proce¬ 
dure to factorize the poles in the regulator e, and the way the subtractions 
are done, is common to both parts. However, contour deformation is only sup¬ 
ported for loop integrals, as the analytic continuation can be performed in 
an automated way, following Feynman’s i6 prescription. More general, user- 
dehned loop integral functions which share the analytic continuation proper¬ 


ties of loop integrals can also be treated in the loop directory, see Section 5.10 
for a detailed example. 


The basic operational sequence of the program is shown in Fig. We also 
show a diagram for the internal SecDec directory structure in Fig. 
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Fig. 1. Flowchart showing the main steps the program performs to produce the 
numerical result as a Laurent series in e. 



Fig. 2. Directory structure of SecDec version 3. 

































































3.2 New features 


3.2.1 New user interface 

The very first change to mention is that the command to lannch SecDec is 
now secdec rather than launch. 


The restrnctnring of the nser interface described below mainly applies to the 
loop setnp. In the general branch(=folder), the program is now also lannched 
by the command secdec, bnt otherwise the nser interface is still the same as 
in version 2. The user interface in the loop branch has been restructured with 
the following aims in mind: 


facilitating scans over large ranges of kinematic parameters 
making the usage of SecDec on a cluster straightforward 
facilitating interfaces to reduction codes like Reduze 
LiteRed 


35 


33 , Fire 34 


or 


• allowing the user to dehne his/her own names for kinematic invariants 

• reducing the mandatory user input to a minimum. 


For these reasons, we have decoupled the dehnition of the numerical values for 
the kinematic invariants from the hie paramloop. input and introduced a sep¬ 
arate input hie kinem. input to dehne the numerical values for the kinematic 
points to be calculated. Each line in kinem. input denotes a new kinematic 
point. The usage is described in detail in Section 4^ and in Appendix A.l[ 


We also changed parts of the Mathematica hie which serves to dehne the 
graph to be calculated. The default name has changed from templateloop .m 
to math.m (however, the user can still give it any name, e.g. mygraph.m). The 
default name for the parameter input hie is now par am. input, the default 
name for the hie dehning the numerical values for the kinematic invariants is 
kinem. input. Accordingly, the command to call SecDec is now either 
secdec 

if the default names for the input hies are used, or 

secdec -p <myparam.input> -m <mygraph.m> -k <mykinem.input>. 
Please note that what has been the “-t” option is now renamed into the “-m” 
option. 


The main changes in the graph dehnition hie math. m are 


The list of propagator powers, powerlist, (also called “indices” in the lit¬ 
erature) can now also contain zero or negative entries, where the latter 
correspond to inverse propagators forming a numerator. Such integrals are 
calculated using the algorithm described in Section 3.2.4 This way the in¬ 
terface to reduction programs providing the master integrals in the form of 










lists of indices is straightforward. 

• The symbols for the kinematic invariants are not predehned, but can be de- 
hned by the user. This also means that the user needs to dehne the expres¬ 
sions for the scalar products of Lorentz vectors (in ScalarProductRules) 
ocurring in the graph. 


The output of SecDec is such that the entire output directory generated by 
the algebraic part of SecDec can be transferred as a standalone archive to 
another machine or cluster, where the numerical evaluation of all kinematic 
points can be submitted in parallel. Optionally, the user can also evaluate 
selected pole coefficients individually, see Section 4^ The output directory 
structure which will be created by SecDec when calculating a loop diagram 
is depicted in Fig. 



■(Makefile kinematics 


.input^ 


^params^ -j^files to extract the parameters defining the integral from the user’s inputJ 


-^file FUN.m containing the graph polynomials U and the numerator'^ 


^decompositionJ - ^contains the results of the decomposition as lists for each pole structureJ 


^subexpj— 


-^files to launch the subtractions and expansions in gj 


[numeric^ - ^produced functions and integration files for each poles structure and each order in cj 


(clusterj— 


-^scripts for job submission to a cluster^ 


^epstothe^ 

(210h0} ;. [e^tothe-1^ 

^epstothe"^ 
[epstothe-1 j 

^epstotheO^ 
( ^OlOho) .[^epstotheO^ 


(results)— 


—(result files, data files for gmiplotj 


(auxres)— 


-(auxiliary files for prefactor, result files for each order in e) 


Fig. 3. Output directory structure generated by the algebraic part of SecDec, 
and example of a pole structure in the numerics folder containing maximally two 
logarithmic poles. 
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3.2.2 Decomposition strategies 

It is well known that the iterated sector decomposition algorithm can rnn into 
an inhnite recnrsion if the variables to be rescaled are chosen in an inconvenient 
way. As an example, consider the fnnction 

f{xi,X 2 ,x^) = xl +xlx^, , ( 6 ) 


and snppose we choose to rescale xi and 0 : 3 . The replacement xi = X 3 ti in the 
snbsector associated with 6 {x^ — xi) leads to f = X3 {X3 1 \ + xf). Snbstitnting 
now X 2 = X 312 in the corresponding snbsector remaps to the original fnnctional 
form, so we generate an inhnite recnrsion. In this simple example we also see 
that the problem is avoided by choosing to rescale Xi and X 2 instead. 


The problem can be stndied systematically by mapping it to a problem in alge¬ 
braic geometry, which in tnrn allows to hnd procednres which are gnaranteed 
to terminate. This was hrst noticed in Ref. by mapping it to Hironaka’s 
polyhedra game, and then rehned in Refs. 


28,36 , see also Ref. 27 


In SecDec-3.0, we have implemented two additional decomposition strategies, 
which are based on the formalism ontlined in Refs, 
guaranteed to stop. 


36 and are therefore 


SecDec- 3.0 uses the program Normaliz 37,38 for the calculation of convex 


polyhedra and their triangulations, needed for the geometric decomposition 
strategies. 


The hrst method G 1 is an implementation of the original algorithm by Kaneko 
and Ueda. Details can be found in Refs. 


36 


The second decomposition algorithm G2, while also based on the ideas of 
Refs. 28,36 , dihers from the original algorithm. Strategy G 2 uses the Cheng- 


Wu theorem ^39 to integrate out the 5-distribution in Eq. ([^, instead of 
the primary sector decomposition employed in other decomposition strategies. 
The Cheng-Wu theorem states that, instead of the 5-distribution in Eq. (f.it 
is also possible to use — ^i)) where 5 is a subset of Feynman parameter 

labels occurring in the integral, provided that the integrations for the remain¬ 
ing Feynman parameters, i.e. those with labels not in 5, are performed from 
zero to inhnity. In our implementation we chose S = {A^}, which amounts to 
setting the parameter xn equal to one, and integrating over xi,..., xat-i from 
zero to inhnity. Even though this breaks possible symmetries among Feynman 
parameters, it turned out to be benehcial with regards to the overall number 
of produced functions. 


Following this procedure, the Minkowski sum of the Newton polytopes of 7/, 
D and the numerator are calculated. The Newton polytope of a polynomial 
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is defined as the convex hull of its exponent vectors. Besides the dehnition as 
the convex hull of a set of points, a polytope A can also be specihed as an 
intersection of half-spaces |40] : 

A = Pi |m G I (m, ni?) -|- > o| . (7) 

F 

Here denotes the primitive normal vector of the facet labeled by F, while 
ap is its distance to the origin. The index F runs over all facets of the polytope 
A. Normaliz implements efficient algorithms to translate between the two 
descriptions. 

For each sector j a variable transformation 

Xi = n (8) 

FeSj 


is performed which remaps the integration variables to the unit hypercube. 
Here the standard basis of is given by the vectors e*. The sets of facets 

Sj are chosen in the following way: 

For each vertex of the polytope, Sj contains the facets incident to it. If the 
vertex lies in more than N—1 facets, the set is decomposed further using a tri¬ 
angulation. The number of generated sectors depends on the way this triangu¬ 
lation is performed. SecDec uses the triangulation algorithm of Normaliz. 


Table [T] compares the number of sectors produced by the decomposition algo¬ 
rithms implemented in SecDec-3.0. The geometric strategy G2 usually gen¬ 
erates the lowest number of sectors. In cases where the integrals have spurious 
poles, which cancel in the hnal result, strategy X can make up for its larger 
number of sectors by producing less spurious poles. 


3.2.3 Improvements in speed 

Towards the goal of using SecDec for the computation of a large number of 
master integrals occurring in amplitudes beyond one loop, it is vital to improve 
on the numerical evaluation times. Several steps towards this goal have been 
undertaken in version 3 of the program. 


CQuad 


As the integrators included in the Cuba library 41 43 are optimized for 


multi-dimensional parameter integrals, it is advantageous to include an ad¬ 
ditional numerical integrator which is dedicated to a very fast evaluation of 
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Diagram 

X 

G1 

G2 

1 1 1 

1 1 1 

282 

266 

166 


1 s 

8s 

4s 


368 

360 

235 


1 s 

9s 

5s 

-r-r ' 




-*-- / 

2 6 '‘'•-.4'' 

548 

506 

304 


3s 

15s 

4s 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

infinite 

72 

76 


recursion 

5s 

1 s 





' 1 -'I'"'' 




r" : r 




'1 1 / 

'v. __ 

27336 

32063 

27137 


5510 s 

11856s 

443 s 

Table 1 


Number of sectors produced by the implemented decomposition strategies and tim¬ 
ings obtained with our implementation of the algorithm. Dashed lines denote mass¬ 
less propagators. 


one-dimensional integrals. One of the currently fastest and still adaptive inte¬ 
grators is Cquad 44 . A version written in the C programming language is 
included in the GSL library 45 . To use it as a standalone tool and to make the 


SecDec package as slim as possible, only the necessary GSL files are included 
in our new release. 


Compilation times 

The setup of the program is such that the algebraic part needs to be performed 
only once for a given diagram, creating functions where the kinematic invari¬ 
ants or other parameters are still in symbolic form. Their numerical values 
only need to be specihed for the numerical integration, without the need to 
redo any algebraic step (except if special kinematics are chosen which would 
change the singularity structure, e.g. massive lines which would lead to IR 
singularities in the limit m —)■ 0 cannot be assigned a zero mass value). 

In previous versions of the program, numerical integration files for each kine¬ 
matic point were generated. This is evaded in the new version, allowing to 
compile the resulting functions once and for all, passing the numerical values 
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as arguments. This has a substantial effect on the performance, as was noticed 
in the inclusion of the program SecDec into FeynHiggs 46 51 to compute 


34 mass configurations of several two-loop two-point integral topologies in an 
automated way, see Ref. 52 . For the latter calculation, a preliminary version 


of the one presented in this paper was used. 


Apart from these changes, we also altered the type of timings displayed in the 
result files. While in previous versions of the program, the processor time spent 
on only the numerical integration part was measured in clock ticks per second, 
we now switched to wall clock times displaying the real time needed for the 
numerical routine (including the presampling when the contour deformation is 
switched on). Even though the real times vary depending on how many other 
processes a user is currently running, we found this information more valuable 
to the user. 


3 . 2.4 Inverse propagators 

As an alternative to the possibility to specify the numerator of a loop inte¬ 
gral as a list of scalar products of loop and external momenta, propagators 
with negative (integer) powers are supported now, which can be convenient 
depending on the way the calculation is organized. 


The implementation of such inverse propagators relies on the following gener¬ 
alization of the Feynman parametrization, which accounts for negative propa¬ 
gator powers by calculating the derivative instead of the integral of the corre¬ 
sponding parameter 10 . For a D-dimensional L-loop integral with N propa¬ 
gators and indices z/i, • • • , z/^r, of which only A),, say the first A),, are positive, 
it reads 


G = (-1)^‘''^F ( W 

N 

n 


,00 


Ui-l 


f)i 




^ 1 - XI 


k=l 







( 9 ) 


^Np + 1 = ---=XN=0 


where Nj, denotes the sum of all positive and negative indices, whereas is 
only the sum of the positive ones. The functions IF and U are defined in the 
usual way, including, however, the propagators with negative powers as well. 
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After performing the derivatives, the result can be written in the form 


G = 


_ ^ _ 

•^0 


LD\ 

,oo 

2 ) 


Af 



-p'i-l 


X, 


Nr, 


ni'i) 


i 1 - It 


k=l 


( 10 ) 


where IAq J-q J^|xjvp+i=---=xjv= 0 ) Nu,n ^j=Np+i 1^*1 

Np,p — Ny, and A/ is a polynomial in a:i, • • • , a^Ar^, which is calculable in terms 
of the partial derivatives of and U up to degree Ny^n. 


The validity of Eq. (§ can be seen easily making a detour over to Schwinger 
parameters. Depending on the sign of its index, one can introduce a Schwinger 
parameter x for a propagator P~^ 


via 


30 


1 

pa 


r(a) ^0 

xP 

dx\°‘\ 


dx x^ 


x=0 


a > 0 
a < 0 . 


( 11 ) 


Performing the usual procedure of completing the square in the exponential, 
shifting and integrating out the loop momenta, it can be shown that the 
Schwinger-parametrized form of G in terms of P and U. is given by 


G=(-l) 


Ny. 


noo 


dx 

i=l 


x: 




N 


r ( D ) ,=^+1 


Q\^j\ 

dx^p ' 


U 


D 

' 2 exp 



3;jVp+i=-"=a;iv=0 

( 12 ) 


This equation can be transformed to Feynman parameters by inserting 


Np 


1 = I dr (5 [ r — 

i=l 


(13) 


and rescaling Xi —)■ rxi for z = 1, ■ ■ ■ ,N afterwards 30 . Finally one performs 
the integral over r using the relation 


roo 

/ dr r'^“^ exp (—rc) = 
Jo 


T{N) 




(14) 


to complete the proof of Eq. §. 


3.2.5 Linear propagators 

As mentioned already, SecDec can also deal with “propagators” which do 
not have a quadratic, but only a linear term in the loop momentum. They 
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occur for example in heavy quark effective theory or in non-covariant gauges. 
An example is given in Section |5.8[ where we calculate 


I 


d^k 1 

ivr? i 5)((fc - piy + i 6){2k ■ V + i 6) 


(15) 


with 7 ^ 0, n ■ Pi = 0. 


Please note that the +i 6 term in the linear propagator {2k ■ v + i6) hxes the 
analytic continuation prescription. The contour deformation implemented in 
SecDec is constructed under the assumption that each propagator carries a 
+i 6 term. 


We should emphasize in this context that great care needs to be taken when 
calculating integrals in the light-cone gauge. In this gauge, one has an addi¬ 
tional linear term q ■ n in the denominator, where q is the loop momentum 
and n is the gauge hxing auxiliary vector with = 0. However, the “naive” 
+i5 prescription is problematic in the light-cone gauge. Therefore the so-called 


Principal Value (PV) piesciipiion 53,54 has been suggested. 


—)■ lim 


+ 


= lim 


q ■ n 


q-n s^o^q-npid q-n — idj s^o {q ■ + 6“^ ^ 


(16) 


but it is not compatible with Wick rotation, which makes the calculations 
based on this prescription quite cumbersome. The only prescription which 
places the poles in the complex q^ plane in the same way as for the prop¬ 
agators in covariant gauges when integrating over the energy component of 
the loop momentum, and therefore preserves the causality behaviour, is the 
Mandelstam-Leibbrandt (ML) prescription 


55-57 


1 1 q-n* 

- —)■ hm - = hm - . 

q-n s^o+q - n + i 6 sign{q - n*) s^o+{q - n){q - n*) + i 6 


( 17 ) 


However, this prescription introduces the dual vector n* {{n*)^ = 0,n-n* ^ 0) 
and therefore leads to a result containing more invariants than in the presence 
of only one auxiliary vector n. 


The above prescriptions are not compatible with the automated contour de¬ 
formation as implemented in SecDec, therefore integrals in the light-cone 
gauge only can be evaluated in the Euclidean region. 


3.2.6 Using Mathematica NIntegrate 

For the numerical integration part, the user of version 3 can choose to either 
use an integrator from the CUBA library [^, or Mathematica’s NIntegrate. 
The available Mathematica options to the latter command can also be specihed 
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directly in param. input. This feature is available in the loop setup, while 
in the general setup, the Monte Carlo integrator Bases 58 is available 


in addition to Cuba. An example for the usage of NIntegrate is given in 
Section 15.91 


3.2.1 Epsilon-dependent dummy functions 

This is a new feature which is available in the general branch. It allows to 
dehne dummy functions which depend on the integration variables, but are not 
explicitly included in the process of iterated sector decomposition. The dummy 
functions are included as external functions in the numerical integration step, 
where the correct replacement of the integration variables, which have been 
rescaled during the iterated decomposition, has been performed automatically. 
This can be very useful in cases where large, but hnite, expressions should be 
integrated together with a function where the poles need to be factorized. 

In version 2 of the program, such a feature was already present, but it did not 
allow to have e-dependent expressions to be masked by the dummy functions. 
Now the dummy functions can depend on e and the program makes sure that 
the various orders in e will be combined automatically with the corresponding 
pole coefficients from the sector decomposition. 


4 Installation and usage 


4.1 Installation 


The program can be downloaded from http://secdec.hepforge.org. Un¬ 
packing the tar archive via tar xzvf SecDec-3.0. tar. gz will create a direc¬ 
tory called SecDec- 3 . 0 . Running make in the SecDec-3.0 directory will call 
the install script, which will check whether Mathematica and perl are present 
and compile the numerical integration libraries Cuba 41 43 , Bases 58 
and Cquad 44 , along with the quasi-random sequence generator Sobol, 


which come with the package. Prerequisites are Mathematica version 7 or 
above, perl (installed by default on most Unix/Linux systems), a C-I--I-03 com¬ 
pliant compiler, and a Fortran compiler if the Fortran part is used. Contour 
deformation is not available in Fortran. 


Please note that the current stable release of Mathematica (vlO.0.2) has a bug 
which causes non-interactive sessions, such as those used by SecDec, to hang 
if parallel kernels are launched. Therefore we have set nbmathsubkrnls=0 
as the default. However, for complicated integrands, using parallelization for 
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the algebraic part can lead to substantial improvements in speed. In such 
cases, we recommend using a different version of Mathematica and enabling 
the parallel execution by setting nbmathsubkrnls to a non-zero number in 
the input parameter file. 


The program NoRMALiz 2.10.1 37, 38p~ is needed for the geometric decompo¬ 
sition strategies G1 and G2. Precompiled executables for different systems can 
be downloaded from http://www.math.uos.de/normaliz/Normaliz2.10.1/ 
and have to be moved into the src/ subdirectory of SecDec-3.0. 


The user can check whether the installation was successful with the command 
make check, which will run a few test examples and compare the results to 
the pre-calculated result coming with the program package. 


4-2 Usage 


In the following we refer to three directory structures: 

• input directory: the directory in which the user’s input files (parameter 
file, math file and kinematics file) reside, 

• output directory: the directory into which SecDec- 3.0 will write output 
files 

• SecDec directory: the location of SecDec-3.0. 

The secdec script is located in the main SecDec directory. It is recommended 
to add the path to the secdec script to the default search paths, so that it 
can be called from anywhere on the user’s system. In the following, we assume 
that secdec was added to the search path, so that it can be called without 
always specifying the path to the script explicitly. 

In a first step, the user should create templates for the input files using the 

command 

secdec -prep. 

The templates will have the default names par am. input, math. m, kinem. input. 
The user should edit these files to define the integral to be calculated and the 
parameters for the numerical integration. If the default names for the input 
files are kept, the program is called from the input directory with the com¬ 
mand 
secdec. 

If the user has renamed the input files, the command is 

secdec -p <parameter_file> -m <math_file> -k <kinematics_file>. 

^ We should note that version 2.10.1 for our purposes works better than the cur¬ 
rently most recent version 2.12.2. 
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The output directory will be created by SecDec. By default the output di¬ 
rectory will be a subdirectory of the input directory and will have the name 
of the graph specihed in the parameter hie. 


An additional option, which can be very useful if the user would like to fac¬ 
torize the poles from functions which do not have the standard Feynman 
parameter integral format including an overall 5(1 — is the user-defined 

option, which is invoked by adding a “-u” to the program call, i.e. 
secdec -u. 

This option skips the primary decomposition step, as it assumes that there is 
no 5-constraint between the Feynman parameters. It requires a certain input 
format for the functions to be decomposed, given in the hie math. m, which is 
explained in the example in Section 5.10 Contour deformation is available for 


these functions, where it is assumed that the user dehned function which takes 
the place of the graph polynomial can be integrated on a contour following 
Feynman’s -\-i6 prescription for the deformation of the contour. Templates for 
the input hies in the user-defined setup can be generated using the command 
secdec -prep -u. 


Likewise, secdec -prep -g will generate templates for the general setup. 
A commented list of all options can be generated by secdec -help. 


4-3 Description of the input files 


par am.input 

The only mandatory helds to be edited by the user in param. input are 
graphname: a name for the integral to be calculated, and epsord: the or¬ 
der in epsilon the user wishes to expand the result. All other parameters take 
default values if not specihed. If a template for param. input is generated 
by secdec -prep, the default values for the options will be hlled in already. 
However, they need not be listed in param. input. The examples in Section]^ 
contain some short input hies where only the helds relevant to the calculation 
of this particular diagram are specihed in the input. All the input options are 
described in detail in Appendix |A.1.5 


math.m 

The hie math. m should contain the dehnition of the graph to be calculated (in 
Mathematica syntax). The template generated by secdec -prep contains the 
dehnitions for a one-loop box diagram. All the settings are explained in detail 
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in Appendix |A.1[ 


kinem. input 

This file is needed for the numerical evaluation to assign numerical values to 
the invariants appearing as symbolic parameters in the functions generated 
during the algebraic part. The symbols for the kinematic invariants have been 
defined by the user in math.m in the lists Kinematicinvariants and Masses. 
Each line in kinem. input should correspond to one kinematic point, where 
the order of the numerical values should match the order of the invariants 
given in these lists, and the masses should always come after the invariants 
formed by Lorentz vectors. Further, each line should start with a label for the 
numerical point (“pointname”), such that the job submission hies and result 
hies for diherent kinematic points belonging to the same graph can be distin¬ 
guished. 

Example: 

Kinematicinvariants = {s,t} 

Masses = {mlsq,m2sq} 

If the user would like to calculate three numerical points, where the Mandel¬ 
stam invariant s can take the values 4,6,8, the value for t is hxed to -0.3, the 
(squared) masses are 1.4 and 2.5, and the points are just labeled by pi ,p2,p3, 
then kinem. input should look like 
pi 4. -0.3 1.4 2.5 

p2 6. -0.3 1.4 2.5 

p3 8. -0.3 1.4 2.5 

Assuming a case in which the user already has tables of numbers for the 

invariants, which do not contain a “pointname” in the hrst row, the script 
addpointname can be used to insert these labels into the data hies. The script 
is located in the folder loop/src/shell. The loop over all kinematic points 
listed in kinem. input is performed automatically. Lines which are commented 
out by a are ignored. 


4-4 Run modes of the program 


The program consists of an algebraic part; factorization of the poles, a nu¬ 
merical part: integration of the pole coefficients, and a collection part: 
collection of the numeric results. If the user only would like to use the alge¬ 
braic part, it is sufficient to issue the command 
secdec -algebraic . 

Similarly, to perform just the numerical integrations, the command 

secdec -numerics 

should be used. The command 

secdec -collectresults 

hnally collects the results. (Any unique substring of the command hags, such 
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as e.g. secdec -c [ollectresults] has the same effect.) 

In single machine mode, all these steps will be done antomatically in one go 
if the user sets the exeflag in param. input to its maximal value (3), which 
is also the default. 

In cluster mode (clusterflag=l), it is not possible to collect the results au¬ 
tomatically, as the information whether all jobs are done is not automatically 
available. Further, the idea is that the user performs the algebraic step once 
and for all on a single machine (where Mathematica is available), and then 
transfers the graph directory to a cluster, where the numerics can be run. 
Therefore, in cluster mode, the three commands secdec -algebraic, secdec 
-numerics, secdec -collectresults have to be issued separately. 

The various possibilities for the user to control the different stages of the 
calculation are shown in Table The optional detailed or basic commands 
must be run in the order presented. If the user specihes an exe flag and calls 
SecDec without a basic or detailed command then all tasks with a lower 
exe flag will also be executed. If the user calls SecDec with a basic or de¬ 
tailed command, for example secdec -subexpand then only that task will be 
performed. Specihcally, in this example makeparams, makeFU, decompose and 
prepare subexpand would not be executed. 


4-5 Evaluation of selected pole coefficients 


In version 3, the user has the possibility to calculate certain pole structures/ep¬ 
silon orders individually. The pole structures are labeled by a string of the type 
i / j hk, where I stands for linear poles, h stands for higher than linear poles, 
while the number i of “usual” logarithmic poles is put at the beginning of 
the string. For example, the pole structure 211h0 means “2 logarithmic poles, 
1 linear pole, 0 higher than linear poles”. By “linear pole” we mean that a 
factorized Feynman parameter has an exponent of the form For more 

details we refer to Ref. (^. 

Each pole structure contains several orders in the e-expansion, ranging from 
the maximal pole of the given pole structure to the expansion order in e speci- 
hed by the user. The coefficients of a certain order in e are stored in subfolders 
of the pole structure folders, which are labeled epstothe-2, epstothe-1, etc. 
The user can select different pole structures as well as individual e orders by 
commands of the form 

secdec [-subexpand/preparenumerics/numerics] [-polestructs= 
<polestructs>] [-epsords=<epsords>]. 


20 


exe- command 
flag (detailed) 


command 

(basic) 


functionality 


> 0 ^ 


makeparams 

makeFU 

decompose 

preparesubexpand 


> algebraic 


subexpand 


> 1 ^ 


preparenumerics 


> 2 


numerics 


numerics 


3 


collectresults 


collect 


extracts parameters from the graph 
definition given in math. m 
constructs the graph polynomials T and U 
performs the iterated sector decomposition 
writes the Mathematica files which are 
needed to run the subtractions and 
e-expansions for each pole structure 
performs the subtractions and e-expansions 
for each pole structure and writes the 
functions f*.cc (resp. f*.m) to be evaluated 
numerically for each pole structure 
writes the files needed to perform the 
numerical integration of each pole structure 
(including the scripts for job submission 
to a cluster) 

performs the compilation and runs 
the executables 

performs the collection of the results 


Table 2 

The different execution stages of SecDec and the possibilities for the user to steer 
them. The optional detailed or basic commands can be issued with, e.g. sendee 
-algebraic, and must be run in the order presented. 

The behaviour of the program is the following: 


• if polestructs is not specified, loop over all (contributing) pole structures 

• if epsords is not specihed, loop over all e orders 

• if polestructs is specihed either as a list 

secdec [-subexpand/preparenumerics/numerics] -polestructs= 
210 h0,110h0, or separately 

secdec [-<subexpand/preparenumerics/numerics>] 

-polestructs=210h0 -polestructs=110h0 

loop over just these pole structures (for all e orders) 
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• if epsords is specified either as a list or separately, loop over just these e 
orders (for all contributing pole structures) 

• if both polestructs and epsords is specihed, loop over just these e orders 
and pole structures. 

If the togetherflag is set to one, all pole structures have been combined into 
one folder called together. In this case, the individual e orders can still be 
calculated separately with the same logic as above, except that there is only 
one pole structure called “together”. For example, in order to perform the 
numerical evaluation only for the 6““^ part of all orders in e contained in the 
together folder, the command would be 
secdec -numerics -epsords=-4. 


5 Examples 


The folders loop/demos and general/demos contain several examples, where 
some of them are tailored to demonstrate a certain new feature, like the us¬ 
age of a numerical integrator from Mathematica instead of the Cuba library 
(NIntegrate), the calculation of tensor integrals (box_withnumerator_2L), 
integrals with e-dependent propagator powers (epsprops_triangle_3L), or 
integrals where strategy X does not stop (geostrategy_2L). 

The input hies for the demos are kept minimal, so the user can see which of 
the input parameters are the relevant ones, while the unspecihed parameters 
take default values. 

The results the user should obtain if the example has run correctly are given 
as hies with the extension . res in the subdirectory check of the corresponding 
example. 

The timings we give for the examples below have been obtained on an Intel(R) 
i7-4790K @ 4.00GHz machine. We give the overall time, i.e. not only the time 
for compilation and numerical integration, but also including the time taken 
by the algebraic part. 


5.1 One-loop box 


This demo is contained in the folder l_box_lL. It contains a simple example 
of a one-loop box graph with one oh-shell leg and one internal mass, see 
Fig. 1^ Two kinematic points are calculated, corresponding to the two lines in 
kinem.input. 
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Pi 


V_ 


P3 


Fig. 4. The one-loop box example, with pi being an off-shell leg. Massless propaga¬ 
tors and light-like legs are shown as dashed lines. 

The overall time to obtain results for the two kinematic points calculated in 
this example was 41 seconds. 


5.2 Two-loop triangle 


The example 2_triangle_2L (denoted by P126 in version 2) is a two-loop 
three-point function containing a massive triangle loop, see Fig. Analytical 
results for this diagram can be found e.g. in Refs. 60,61 , a threshold scan is 
given in Ref. [^. 



Fig. 5. The graph P126, containing a massive triangle loop. 
The overall time taken by this example was 144 seconds. 


5.3 Two-loop non-planar box with internal masses 


This example is contained in the folder 3_nonplanarbox_2L. 


It is a 7-propagator non-planar two-loop box integral where all propagators 
are massive, using mi = m 2 = m^ = m^ = m, m 3 = m 4 = my = M, 
p\ — pi — p\ — The labelling is as shown in Fig.[^ 


Numerical results for this integral have hrst been calculated in Ref. 62 using 


a method based on extrapolation in the i5 parameter. We give results for 
m = 50, M = 90,S23 = —10“^. A scan over the invariant S 12 can be found in 
Ref. 


The overall time taken by this example was 186 seconds. 
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Fig. 6. Labeling of the masses for the non-planar graph . 
Two-loop planar box with loop momenta in the numerator 


Pi 


P2 


ki 


ki - k2 


Pi 


P3 


Fig. 7. Momentum labelling for the two-loop planar box. All propagators and ex¬ 
ternal legs are assumed to be massless. 

In this example, which can be found in the folder 4_box_withnumerator_2L, 
we consider a planar massless two-loop box diagram with an additional prop¬ 
agator that is raised to a negative power, i.e. is in the numerator. 

Using the momentum labelling shown in Fig. we choose 

ih + psf = kl + 2ki ■ ps (18) 

as the inverse propagator. There are two ways to compute the integral with 
this numerator: 

• Add {ki to the propagator list and specify its index as —1 in the 

power list. This method is demonstrated in the default input files in the 
example folder and can be run by just typing secdec. 

• Calculate two integrals, with kf and 2ki ■ ps, respectively, as numerators. 
In the first integral, the numerator can be cancelled against a propagator. 
This pinched graph can be calculated by simply removing this propagator 
from proplist, or by setting its index to 0 in the powerlist. For the second 
integral, we explicitly specify the numerator in the mathfile math.tensor .m. 

These integrals can be run by specifying the corresponding input files: 
secdec -m math_pinched.m -p param_pinched, input 
secdec -m math.tensor .m -p param.tensor. input 

We emphasize that the two methods to calculate integrals with numerators 
cannot be combined, i.e. when an index smaller than zero is specified, the 
numerator given in math.m must not include any loop momenta. 
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The overall time taken by the example was 14 minutes for the graph with 
inverse propagator, i.e. the default input hies param. input, math.m, below 3 
minutes for the graph with specihed numerator and 19 seconds for the pinched 
graph. 


5.5 Two-loop pentagon 


This example is contained in the folder 5_pentagon_2L. The calculation of 



\ 

\ 


^ 4 

^- 



Fig. 8. The numbers label the vertices of the massless two-loop five-point integral. 

the diagram in Fig. [^demonstrates the applicability of SecDec to examples 
with more than four external legs. Assuming the external momentum pi to be 
situated at vertex f, and choosing the invariants formed by two adjacent legs 
as a basis, the following scalar products need to be dehned 

ScalarProductRules = { 

SP[pl,pl]^0, SP[p2,p2]^0, SP[p3,p3]^0, SP[p4,p4]^0, 
SP[p5,p5]^0, SP[pl,p2]^sl2/2, SP [pi ,p3] ^ (s45-sl2-s23)/2, 

SP[pi,p4] (s23-s51-s45)/2, SP[pi,p5] ->s51/2, 

SP [p2,p3] ^s23/2, SP [p2 ,p4] ^ (-s23-s34+s51)/2, 

SP [p2 ,p5] ^ (s34-sl2-s51) /2, SP [p3 ,p4] ^s34/2, 

SP[p3,p5]—j-(sl2-s34-s45)/2, SP [p4,p5] —>s45/2 }; 

which are formed from relations among the kinematic invariants Sij = {pi -|- 
The overall time taken for this example was 9 minutes and 43 seconds. 


5.6 Two-loop with geometric decomposition strategy 


This example is contained in the folder 6_geostrategy_2L. It calculates the 
two-loop box diagram shown in Fig. using the strategy G2. With strategy 
X the decomposition runs into an inhnite recursion. Numerical results for this 
diagram in the Euclidean region have hrst been given in Ref. [^. Analytical 
results can be found in Ref. 


64 . Here we expand the result up to order 
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Fig. 9. The two-loop box example with four massive on-shell legs. Dashed lines 
denote massless propagators. 

The overall time taken for this example was 41 seconds. 


Please note that Normaliz 37,38 must be installed to run this example. 


5.7 Three-loop triangle 


This example is contained in the folder 7_epsprops_triangle_3L. Apart from 
being a 3-loop example, it also demonstrates the usage of e-dependent prop¬ 
agator powers. It calculates the diagram shown in Fig. IT, where one of the 
propagators is raised to the power 1 -|- e. We expand the result up to order e^. 


< 



Fig. 10. The three-loop vertex diagram Agp with the dotted propagator raised to 
the power 1 -|- e. 


The analytical result for this diagram with general propagator powers is given 
in Ref. 65 and is also given in the hie AGlanalytic .m to allow comparisons 


between analytical and numerical results for arbitrary propagator powers. 


The overall time taken for this example was 13 seconds. 


5.8 Linear propagators 


This example is contained in the folder 8_linearprop_2L. It contains the in¬ 
tegral 


I 


d^k 1 

ivr? - pi)2 + i 6){2k ■ V + i 6) 


( 19 ) 


with 7 ^ 0, n ■ Pi = 0. 
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Please note that in the presence of linear propagators, proplist must be given 
in the form of explicit momentum flow, i.e. constructing the graph polynomials 
W from labeled vertices is not possible. 

The overall time taken for this example was 30 seconds. 


5.9 Using NIntegrate instead of Cuba 


This example is contained in the folder 9_NIntegrate_box_lL. It is identical 
to the one-loop box described in Section 5.1, except that the NIntegrate 
function from Mathematica is used for the numerical integration. 


The overall time taken for this example was 40 seconds. 


5.10 User defined integrands 


This example is contained in the folder 10_userdef ined_triangle_lL. In this. 



•" 3 


Fig. 11. The numbers label the vertices of the one-loop three-point integral. Dashed 
lines denote massless internal or external particles, full lines denote massive ones. 

a one-loop three-point function with one massive propagator and one off-shell 


leg, compare Fig. 11, is computed using the user-dehned setup. 


The user-dehned setup supports functions of the type 


Guser = / n , 

0 J = 1 


( 20 ) 


where n is the number of integration parameters Zj to be integrated from 
zero to one and A/" is a numerator function which cannot contain singularities. 
lA and 5F may contain singularities and, if needed, sector decomposition can 
be applied to factorize these. Those parts of a user function which need an 
analytical continuation a la loop integral, should be written to the 5F part, as 
the integration contour is deformed based on the properties of this function. 
It should be noted that, comparing with Eq. (|^, Eq. (20) does not contain a 
^-distribution. 
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In the standard Feynman loop integral case, the 5-distribntion is integrated 
ont dnring the primary sector decomposition. Performing the primary sector 
decomposition of the one-loop triangle diagram, see Fig. 11, by hand, the 
resulting functions read 


Gl,user = J dZidZ 2 Z^^~^{l + Z2)~'^^^^{-SZ2 + m^(l 2:1 + Z2))~^~^ , 

0 

(21a) 

1 

G2,usei = J d^id^2^2 ^ ^(1 + + 2 : 2 ) s-|-m^(l-f ^1-|-Z 2 )) ^ (21b) 

0 

1 

Gs.user = y d^idz2(l + + ^ 2 ) —S^2 +’^^(l + + ^ 2 )) ^ ■ (21c) 

0 


It must be noted that the standard loop prefactor — LD/2) in 

Eq. ([^ is not included in the setup for the userdehned functions. In order to 
get the same result as in the Feynman loop integral case, one therefore has to 
compute 


Gioop — — F(1 -|- e)(Gi ,user “1“ 1-^2,user “1“ fj3,user) 


In SecDec, the functions of Eqs. (21) are dehned in the functionlist 
functionlist = { 

{1, {-l-eps,0}, 


{{(U[z]/.z[l]-^l)/.z[3]->t[l] , expou, A}, 
{-s*z[2] + msq*(1+z[1]+z[2]), expof. A}}, Num}, 
{2, {0,-1-eps}, 


{{(U[z] /. z [2] —>1) / .z [3] —7-t [2] , expou, A}, {-s + msq* (1+z [1]+z [2]), 


expof. A}}, Num}, 

{3, {0,0}, 

{{U[z] /. z [3] —7-1, expou. A}, {F [z] /. z [3] —>1, expof, B}}, Num}}; 


where the functions F [z] and U [z] corresponding to the two graph poly¬ 
nomials of the Feynman one-loop triangle integral, their exponents expou, 
expof and the numerator Num of the integral are separately written to the 
math. m hie using Mathematica syntax 
Num = 1; 

expou = -l+2*eps; 
expof = -1-eps; 

U[z_] := z[l] + z[2] + z[3] 

F[z_] := z [3] * (-s*z [2] + msq*(z[l] + z [2] + z [3])) . 

The user functions are labeled by their hrst entry in the functionlist. If two 





functions share the same exponents expou and expof, and are labeled with the 
same number, they are summarized and then decomposed as if they belonged 
to one sector. 

The second entry of each user function in the functionlist is a list of exponents. 
Here, the exponents of each integration parameter should be given. The hrst 
function Gi^uger of this triangle example contains a factor but no factor 

of Z 2 - This is translated into the statement {-l-eps,0}. 

The third entry of each user function in the functionlist contains a list of infor¬ 
mation on the functions corresponding to the hrst and second graph polyno¬ 
mial after primary sector decomposition in our example. This is the function 
itself, its exponent and the hag A or B. While A denotes that no further decom¬ 
position is needed, B leads to an iterated sector decomposition of the fnnction. 

While the inpnt hies param. input and kinem. input are the same for the 
Feynman loop integral and the userdehned setup, the math.m input hie has to 
contain the above mentioned functionlist dehnition. Unlike the loop setnp, def¬ 
initions for momlist, proplist, numerator, powerlist, ExternalMomenta, 
externallegs, ScalarProductRules need not be inclnded if these are not 
necessary for the dehnition of the integrand in functionlist. 

To run this example, go to the demo directory 10_userdef ined_triangle_lL 
and issne the command secdec -u, the -u activates the userdehned setnp. 

The overall time taken for this example was 7 seconds. 


5.11 Dummy functions 


This example demonstrates the usage of e-dependent dummy functions in 
the calculation of general parametric integrals. It is located in the folder 
general/demos. The integrand of this example reads 

f = {zi + Z 2 )~^~^" dnm}+"( 2 :i, 2 : 2 , 2 : 3 , 2 : 4 ; e) (l\iml~^\z 2 , z^) cnt( 2 : 3 ), ( 22 ) 

with 

dnmi(a, b, c, d; e) = a'^ + + c'^ + + 4:abcd + 2 — ofb^c^df’ + 5ae, (23a) 

dum 2 (a, b) = of + b'^ -\- 0^ + Aab — \J abf -f 3a^6^, f = 0.5 (23b) 

cnt ( 2 : 3 ) = 0 ( 2:3 — C), C = 0.1. (23c) 

Note that the dnmmy fnnctions must not lead to additional singnlarities. 
Before secdec can be lannched to calculate the integral Jq d 2 :id 2 : 2 d 2 ; 3 d ^4 /, 
Fortran code for the calculation of the dummy functions has to be generated. 
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This can be done by specifying the names of the dummy functions, as well as 
their highest power in e, in the parameter hie and running 
perl createdummyfortran.pl -d demos/ -p dummyfunction.input 
from the general folder. This will translate the Mathematica hies contain¬ 
ing the dehnitions of the dummy functions to Fortran code. However, the hie 
cut. f has to be edited afterwards to account for the 0-function. These mod- 
ihcations can be found in the hie cut_edited.f and should be copied to the 
hie cut. f. 

The integration is launched with 

secdec -g -p dummyfunction.input -m dummyfunction.m 

and will take into account the e-dependence of the dummy functions. When 

doing the algebraic and numeric part of SecDec separately, it is possible 

to modify the dummy functions after the algebraic part as long as their e- 

dependence does not exceed the maximum power specihed in the parameter 

hie. Further details can be found in a readme hie located in the general/demos 

directory. 

The overall time taken for this example was 3 seconds. 


6 Conclusions 


We have presented version 3.0 of the program SecDec, which is publicly 
available at http://secdec.hepforge.org, 

The part of the program which allows to calculate multi-loop integrals for 
arbitrary kinematics has been improved in various respects: it contains two 
additional decomposition strategies which are guaranteed to stop, based on a 
geometric algorithm. In addition, it can deal with tensor integrals in the form 
of inverse propagators, or, more generally, can take lists of indices where the 
indices can also be negative or zero. Integrals containing linear propagators 
can also be calculated. The timings also have been improved, and the possibil¬ 
ity to do the numerical integrations for a large range of kinematic points on a 
cluster has been made much more straightforward. Further, it is now possible 
to use an integrator from Mathematica in addition to the Cuba library. To¬ 
gether with the re-structuring of the code, the whole interface has been made 
more user-friendly. 

The part of the program which allows to factorize poles from parameter in¬ 
tegrals which are not related to loop integrals also contains a new feature, 
allowing to introduce some dummy functions which will not be decomposed, 
but can themselves depend on e. 

With all its new features and the new user interface, we believe that SecDec- 
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3.0 will certainly be useful for a number of calculations beyond NLO, in par¬ 
ticular in cases where several mass scales are involved, and where analytical 
calculations of the loop integrals are at their limit. 
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A User Manual 


SecDec has various setups: one for standard loop integrals, termed ‘loop 
setup’ in the following, and one for general parametric functions, termed ‘gen¬ 
eral setup’, as reflected by the two directories loop and general. While the 
parametric functions treated in the general folder are not accessible to con¬ 
tour deformation, the latter is still available for integrals which go beyond 
standard loop integrals. Such non-standard integrals can be dehned by the 
user, and therefore this setup is referred to as ‘user-dehned setup’, which is 
an option within the loop directory. 

The program is called by invoking the script secdec, located in the main 
SecDec directory. We recommend to add the path to the secdec script to 
the user’s default search paths, so that it can be called from anywhere on 
the system. In the following, we assume that secdec was added to the search 
path, so that it can be called without always specifying the path to the script 
explicitly. 


A.l Overview of usage in the loop setup 
A. 1.1 Basic usage 

• Greate templates for the input hies with the command 
secdec -prep 

This command generates the hies par am. input, kinem. input and math.m. 
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• Input the graph name and desired order in e (epsord) into par am. input, 
dehne the loop integral in math. m and specify one or more kinematic points 
in kinem. input. Please note that the order of the values for the kinematic 
invariants given in kinem. input must match the order of the invariants 
given in math.m, i.e. the entries of the lists Kinematic Invariants, Masses 
(masses always listed last). 

• To run SecDec: simply issue the command 
secdec 

If the input hies have been renamed, the command is 

secdec -p <myparam.input> -m <mymath.m> -k <mykinem.input>. 

• Switch to the output directory which has been created. It carries the graph 
name specihed in param. input. The results are in the results folder, con¬ 
taining a hie graph.pointname. res for each point specihed in the kinemat¬ 
ics hie, and plotf ile<i> .gpdat hies for each e-order i, where the results 
for all kinematic points have been appended. 


A. 1.2 Intermediate usage 

The diherent stages of SecDec can also be run separately. Creating and 
editing the input hies as before, one can run the algebraic part only by 
secdec -algebraic. 

This allows, for example, to get a hrst idea of the pole structure generated 
during the decomposition. One can also examine the IF and lA functions in the 
folder named FU in the output directory. 

The numerical part is run by the command 
secdec -numerics. 

The results are collected by 
secdec -collectresults. 

Options: 

• to add a few more kinematic points, for example in the threshold region, 
one can edit the kinem. input hie in the input directory. The user should 
delete the old (already calculated) points from the kinem. input hie and 
add the new points with new point names. 

• One can re-run the numerics in cluster mode by uploading the output di¬ 
rectory to the head node and running the submission scripts in the cluster 
folder. (Note: please examine the submission scripts before submitting them 
to your cluster, they may need to be adapted to the particular cluster 
setup). To collect the results, one should copy the output directory back 
to the original machine (for example by using rsync) and then run: secdec 
-collectresults. 
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A. 1.3 Expert usage 

The various stages of SecDec can be subdivided even more and can be run 
separately by issuing the detailed commands listed in Table For example, 
to only extract the graph parameters and construct the graph polynomials IF 
and W, one has to issue the commands 
secdec -makeparams and then 
secdec -makeFU. 


The user can select pole structures and e orders to be computed. For example, 
the command: 

secdec -numerics -polestructs=210h0,llOhO -epsords=-2,-l 
would compute only the contribution of the 210h0 and llOhO pole structures 
to the and e~^ poles. 

For all possible commands we refer to Table Please note: if an exe flag is 
specihed and secdec is called without a basic or detailed command then all 
tasks with a lower exe flag will also be executed. If the user calls secdec with 
a basic or detailed command, for example secdec -subexpand then only that 
task will be performed. 


A. 1.4 Description of the input files 

The description below is for loop diagrams; the input hies in the subdirectory 
general to compute more general parametric functions did not change be¬ 
tween versions 2 and 3, and therefore we refer to descriptions in Refs. |2 
for the general branch. 

• param. input: (text hie) 

The mandatory parameters the user needs to specify are 

■ graphname: a name for the graph to be calculated 

■ epsord: the desired expansion order in e. 

All other values take defaults if not specihed. A detailed description of all 
options is given is Section |A.1.5[ 

• math. m: (Mathematica syntax) 

This hie contains the dehnition of the graph to be calculated, 
momlist: list of loop momenta 

proplist: list of propagators, either in momentum how representation or 
as a list containing the propagator mass and the labels of the vertices the 
propagator is connecting. Please note that the label of the vertex which 
contains the external momentum pi should be i. For vertices involving only 
internal lines, the labelling is arbitrary. For more details we refer to Ref. . 
numerator: there are two possiblities to dehne a numerator: 

(a) give a list of loop momenta contracted with either external momenta or 
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loop momenta. Each Lorentz contraction should be denoted by while 
each contracted factor should form an element of a list. E.g., for 2 ki-pi k 2 -p 2 , 
the syntax is numerator=2,kl*pl,k2*p2. 

(b) dehne an additional propagator and specify a negative index in powerlist 
(see below). 

The default is numerator=l (scalar integral). 

powerlist: list of porpagator powers, also called “indices” in the literature. 
Can also take zero or negative integer values. 

Dim: the dimension of the loop momenta. The default is Dim=4 — 2e. 
prefactor: the prefactor specihed here will be factored out of the numerical 
result. This means that the numerical result will be divided by the prefactor 
given here. Please note that a factor — LD/2), which comes 

from the Feynman parametrisation, will be included by default in the nu¬ 
merical result, according to the integral dehnition in Eq. ([^. 
ExternalMomenta: list of external momenta occurring in the graph dehni¬ 
tion (only necessary if the graph dehnition contains the momenta explicitly, 
i.e. in momentum how representation). If the length of this list is diherent 
from the number of external legs, give also the number of external legs as 
externallegs=.... 

Kinematicinvariants: list of symbols for kinematic invariants (formed 
from Lorentz vectors) occurring in the diagram. The symbols can be chosen 
by the user. 

Masses: list of symbols for the masses. In the case of massive on-shell lines, 
i.e. p^ = m?, where m is a propagator mass as well as the mass of an external 
leg, the way to proceed is to dehne SP[p,p] —» m2 in ScalarProductRules, 
such that no extra symbol needs to be specihed for p^. 

Note that the list of masses should always contain the propagator masses 
used in proplist. 

ScalarProductRules: list of replacement rules for the kinematic invariants 
formed by external momenta. 

splitlist: allows to specify those Feynman parameters which may have 
an endpoint singularity at Zi = 1, such that the integral will be split at 1/2 
and the singularity at one will be remapped to a singularity at zero. 

• kinem. input: (text hie) 

Should contain numerical values for the kinematic invariants, in the same 
order as the symbols for the kinematic invariants given in the helds 
Kinematicinvariants and Masses in the Mathematica input hie (called 
math.m here). The numerical values for the masses should always be listed 
after the invariants formed from Lorentz vectors. 

If kinemloop.input contains several lines, each line will be evaluated as a 
new kinematic point. In order to be able to distinguish the runs/results for 
the diherent kinematic points, each array of numerical values should have 
a label prepended which dehnes the “pointname”. Example: if mathloop.m 
contains Kinematicinvariants = {s,t,plsq} and Masses={mlsq}, the cor- 
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responding kinematics input file where mlsq takes the value 1 for point pi 
and 5 for point p2 should look like (the values for s,t,plsq are e.g. 500, -88, 
100 ) 

pi 500 -88 100 1 
p2 500 -88 100 5 


A.1.5 Detailed description of all options in param. input 

It should be emphasized again that the only mandatory helds are graph and 
epsord, all other parameters take default values if not specihed, which are 
given in brackets after the keyword in the description below. 


graph The name of the diagram or parametric function to be computed is 
specihed here. The graph name can contain underscores and numbers, but 
should not contain commas. 

epsord The order to which the Laurent series in e should be expanded, start¬ 
ing from The value of epsord=0 will calculate the pole coefficients 

and the hnite part. Note that epsord can be negative if only the pole coef- 
hcients up to a certain order should be computed, 
outputdir 0 specihes the name of the directory where the produced hies 
will be written to, the absolute path should be given. If left empty, a 
subdirectory of the input directory with the name of the graph will be 
created. 

contourdef (False) The contour deformation can be switched on or oh by 
choosing contourdef= True/False (lower case letters are also possible). For 
multi-scale problems, respectively diagrams with non-Euclidean kinematics, 
set contourdef = True. 

lambda (1.0) A is a parameter controlling the contour deformation. The pro¬ 
gram takes the A value given by the user in the input hie as a starting point. 
The program then performs checks and optimizations to hud an ‘optimal’ 
value for A. The user should pick an initial value which is rather large, as the 
program will decrease A appropriately, while it cannot increase A. Values of 
A between 1 and 3 are usually a good choice, 
strategy (X) Choice of the decomposition strategy. The default is X, which is 
the same strategy as in previous versions, which is usually the most efficient 
one, but is not guaranteed to stop. If the decomposition does not seem to 
stop, strategies G1 or G2 should be chosen, which are based on a geometrical 
algorithm as described in section 3.2.2, and which cannot run into an inhnite 
recursion. Within strategy G2, no primary sector decomposition is done. 
Therefore it is obviously not possible to specify only selected primary sectors 
to be calculated when using this strategy, 
exeflag (3) The exeflag can be used to execute the program only up to a 
certain stage. There are three basic stages: i) the algebraic part, ii) the 
numerical part and hi) collecting the results. The algebraic part can be 
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further split into a part doing only the iterated sector decomposition (e.g. to 
get an idea about the pole structure), and a part performing the subtractions 
and the expansion in epsilon. The values of the exeflag correspond to the 
following stages: 

• 0: The parameters like the number of loops, propagators etc are extracted 
from the user’s Mathematica input hie; the graph polynomials T and hi 
(and the numerator in the case of tensor integrals) are constructed; the 
iterated sector decomposition is done; the scripts subandexpand*. m in 
the graph subdirectory for the subtractions and epsilon expansions are 
created, but not run. 

• 1: The subtractions and expansions in e are performed and the resulting 
functions for the pole coefficients are written to C++ or Mathematica 
hies; all the other hies needed for the numerical integration are created as 

well. 

• 2: Compilation and numerical integrations are performed. 

• 3: The results are collected. 

All exehags imply that the steps corresponding to lower exehags will au¬ 
tomatically be performed as well. However, there is also the possibility to 
skip previous steps by calling SecDec with a basic or detailed command, 
as explained in Section [4)4 For example, to run only the numerical part the 
call is secdec -numerics. This command will not restart the whole alge¬ 
braic part, but just compile the functions and run the executables. An error 
will be thrown if the functions have not been produced beforehand. Table 
gives a schematic overview of the various options to control the program 
how. 

Please note that if the clusterf lag is switched on, it is assumed that the 
user will produce the functions locally and then compile and run them on a 
cluster. Therefore, in cluster mode, the program will perform the algebraic 
part and produce the submission scripts, and then stop, independent of the 
value of the exeflag, as the user should control the job submission (secdec 
-numerics) and the collection of the results (secdec -collect) in cluster 
mode. 


Advanced usage 

togetherflag (0) This hag dehnes whether to integrate subset^of functions 
contributing to a certain e-order separately {togetherflag=0), or to sum all 
functions for a certain order in e prior to integration {tog ether flag=1). The 
latter will avoid large cancellations between results for functions descending 
from diherent pole structures and thus give a more realistic error. However, 
togetherflag=1 is not recommended for cases where the individual functions 

^ The subsets are naturally formed by the fact that functions contributing to a 
certain e-order can descend from different pole structures. 
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are already very complicated. 

grouping (2000) It can be beneficial to first snm a few fnnctions before in¬ 
tegrating them. Choosing a valne for the gronping which is nonzero defines 
an npper limit (in kilobytes) for a file containing the snm of a nnmber of 
fnnctions. The nnmber of kbytes is set by grouping=#khytes. Setting group- 
ing=0, all fnnctions f * . cc resp. f* .m are integrated separately. In practice, 
gronping =0 has proven to lead to faster convergence and more accnrate re- 
snlts in most cases. However, for integrals where large cancellations among 
the different fnnctions occnr, the gronping valne shonld be chosen 7 ^ 0. The 
log files ^results* . log in the resnlts directory contain the resnlts from the 
individnal snb-sector integrations. These files are nsefnl to spot cancellations 
between the individnal fnnctions and to adjnst the settings accordingly. 
IBPflag (0) Set IBPflag=0 if the integration by parts option shonld not be 
used and IBPflag=l if it should be used. IBPflag=2 is designed to use IBP 
relations when it is deemed efficient to do so. Using the integrations by 
parts method takes more time in the subtraction and expansion step and 
generally results in more functions for numerical integration. Its usefulness 
is mainly for cases where (spurious) linear poles of the type are found 

in the decomposition, as it reduces the power of x in the denominator, 
infinitesectors () This field should be empty if the default strategy X or 
the strategies G1 or G2 are applied. It offers the possibility to use a different 
‘heuristic’ strategy for certain primary sectors if they seem to suffer from 
infinite recursion. The primary sectors given in the list (comma separated) 
will then be decomposed with this heuristic strategy. It can avoid infinite 
recursion in cases where strategy X fails, but is not guaranteed to stop. 
For example, infinitesectors=2,3 results in the application of this heuristic 
strategy to primary sectors 2 and 3. In examples with massive propagators, 
one should put the labels belonging to the massive propagators into the list, 
primarysectors () This field allows to calculate selected primary sectors 
only. If left blank, primarysectors defaults to all, i.e. primarysectors=l,...,N 
will be assumed, where N is the number of propagators. This option is useful 
if a diagram has symmetries such that some primary sectors yield the same 
result. It cannot be used in combination with strategy G2. 
multiplicities () Specify the multiplicities of the primary sectors listed above. 
List the multiplicities in same order as the corresponding sectors above. If 
left blank, a default multiplicity of 1 is set for each primary sector, 
rescale (0) If there are large differences in magnitude in the kinematic in¬ 
variants occurring in a diagram, it can be beneficial to rescale all invariants 
by the largest one in order to reach faster convergence during the numer¬ 
ical integration. The rescaling can be switched on with rescale=l (default 
is rescale=0). Please note: If switched on, it is not possible to set explicit 
values (numbers) for any non-zero invariant in the ScalarProductRules= 
conditions in the Mathematica file math. m. 
nbmathsubkrnls (0) Maximal number of Mathematica kernels to be used 
by Mathematica (the iterated decomposition and subtractions are not par- 
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allelized by default). 

smalldefs (0) smalldefs=l minimizes the deformation of the contour. It can 
be useful for example in the presence of oscillatory integrands. 

largedefs (0) If the integrand is expected to have (integrable) endpoint sin¬ 
gularities at Xj = 0 or 1, largedefs=l, can help to have a sufficiently large 
deformation close to the endpoints. The default is largedefs=0. 

Note that setting both flags largedefs and smalldefs to zero is perfectly pos¬ 
sible, as the flags operate on different parts of the deformation internally. 
For more details we refer to Ref. [^. 

optlamevals (4000) The number of pre-samples to determine the optimal 
contour deformation parameter A can be chosen by assigning a number to 
optlamevals. 


Parameters related to the numerical integration and external libraries 

compiler (gcc) Choice of the C-compiler. 

CCargs (-0) Compiler options for the C-compiler to compile the numerical 
integration hies. 

sobolpath 0 The path to sobol can be specihed here, if different from the 
default [path_to_secdec]/src/sobol. The sobol quasi-random number 
generator is only used if eontourdef=True. 

cquadpath () The path to cquad can be specihed here, if diherent from the 
default [path_to_secdec] /src/cquad. Note that the program will use this 
integrator automatically if an integral or a pole contribution is found to 
depend on only 1 Feynman parameter, irrespective of what has been chosen 
as integrator below. 

integrator (3) The program for the numerical integration can be chosen 
here. Vegas {integrator=l), Suave {integrator=2), Divonne {integrator=3) 
and Cuhre {integrator=4) are part of the Cuba library. To choose a numer¬ 
ical integrator included in Mathematica, integrator=5 can be chosen. 


NIntegrate parameters 

NIntegrateOptions (AccuracyGoal->3) Options for the Mathematica NIn¬ 
tegrate command, if integrator=5 is chosen. Example: 
NIntegrateOptions=AccuracyGoal —?• 2,WorkingPrecision—;-12, 

Method—>■ “AdaptiveMonteCarlo” 

Please note that when using NIntegrate, it is not possible to obtain an error 
estimate. 
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Cuba parameters 


Below we describe only a few Cuba parameters, for more information the user 
is referred to the Cuba documentation, Refs. 41,43 . As we have set Divonne 
to be the default integrator, we give some more details for Divonne settings 
below. 


cubapath () The path to the Cuba library can be specihed here, if different 
from the default [path_to_secdec] /src/Cuba-4.1. 
cubacores (1/0) The maximal number of cores Cuba is allowed to use. In 
cluster mode, the default is 1, in single machine mode the default is zero, 
which means that Cuba will use all available idle cores, 
maxeval (10000000) The maximal number of evaluations to be used by 
the numerical integrator. For this held and the helds mineval, epsrel, 
epsabs below, a value can be specihed for each order in e, separated by 
commas and starting with the leading pole. If only one value is given, it 
will be used for all pole coefficients. If the list is shorter than the number of 
orders in e, the last value of the list will be repeated as often as necessary, 
mineval (0) The number of evaluations which should at least be done before 
the numerical integrator returns a result, 
epsrel (l.e-2) The desired relative accuracy for the numerical evaluation, 
epsabs (1. e-6) The desired absolute accuracy for the numerical evaluation. 
These values are particularly important when either the real or the imagi¬ 
nary part of an integral is close to zero, 
cubaflags (2) Sets the Cuba verbosity hags. The default is 2, which means 
that the Cuba input parameters and other useful information, e.g. about 
numerical convergence, are written to the log hie of the numerical integra¬ 
tion. 

seed (0) The seed used to generate random numbers for the numerical inte¬ 
gration with Cuba. The default is seed=0\ Cuba will use the Sobol (quasi-) 
random number generator 

Divonne speeific parameters 

keyl (1500) Determines the number and type of sampling to be used for the 
partitioning phase in Divonne. With a positive keyl diherent from 7,9,11,13, 
a Korobov quasi-random sample of keyl points is used. 
key2 (1) Determines the number and type of sampling to be used for the 
hnal integration phase in Divonne. With 77-2 = \key2\ < 40, the number of 
sampling points is n 2 nneed, where Uneed is the number of points needed to 
reach the prescribed accuracy as estimated by Divonne from the partitioning 
phase. 

keys (1) Sets the strategy for the rehnement phase in Divonne. Setting key3=l 
(default), each subregion is split once more, 
maxpass (4) Number of iterations that are performed before the partition 
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is accepted as final, to pass on to the main integration phase, 
border (l.e-6 ) The points in the interval [0,border] and [1 — border, 1] are 
not included in the integration but are extrapolated from samples in the 
interior. This can be useful if the integrand is known to be peaked close to 
endpoints of some integration variables, 
maxchisq (1.) The maximally allowed a particular subregion is allowed 
to have in the hnal integration phase. Regions which fail this test and whose 
sample averages differ by more than mindeviation move on to the rehne- 
ment phase. 

mindeviation (0.15) Determines if a region which failed the f^st is treated 
further (see maxchisq above). 


Parameters related to the cluster mode 

clusterflag (0) Determines how jobs are submitted. Setting clusterflag=0 
(default) the jobs will run on a single machine, with clusterflag=1 the jobs 
will run on a cluster (the corresponding hies to submit jobs to a cluster will 
be created, see below). 

batchsystem (0) Chooses a format for the scripts steering the submission 
to a cluster. If batchsystem is set to 0, the setup is for the PBS (portable 
batch system). If the hag is set to 1, a user-dehned setup is activated. 
Currently this is the submission via condor, but it can be easily adapted to 
other batch systems by editing the templates in loop/src/numerics/ and 
loop/perlsrc/makej ob.pm. 

clusteroptscompile () In cluster mode: command line arguments passed 
verbatim to the job submission script for compilation jobs on a cluster 

clusteroptsrun () In cluster mode: command line arguments passed verba¬ 
tim to the job submission script for numerical integration jobs on a cluster. 


Parameters related to plotting 

xplot (1) This option can be used to control the format of the data hies where 
the results for a range of kinematic points are listed. The variable xplot 
denotes a position in the list of invariants. The corresponding invariant then 
will be the one which will be plotted on the x-axis. Example: the invariants 
are s,t,u,mlsq,m2sq. If a scan over mlsq has been performed, such that 
mlsq should be plotted on the x-axis, then xplot=4 would tell SecDec 
to write the values for mlsq into the hrst column of the *.gpdat hie. The 
*. gpdat hies produced by the program have the form 
[invariant chosen by xplot] real_result real_error imag_result 
imag_error timing. 

For 3D plots: if xplot is a list of length L, the hrst L columns of *.gpdat 
will correspond to the values of the invariants singled out by the xplot labels 
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(e.g. xplot=l,2 would produce data files for a 3D plot in s,t)- 


A. 1.6 User-defined setup 


This setup allows to define functions in the Mathematica input file 
(math_userdefined.m) which are not standard loop integrals. It is invoked by 
the option -u when calling secdec. While the param. input file has the same 
form as for standard loop integrals, the hie math [userdef ined] .m should 
contain the specihcation of the user-dehned functions. Most helds are the 
same as in math [loop] .m. However, instead of the dehnition of a graph, the 
user can dehne a list of functions to be decomposed. The detailed format is 
specihed in the example described in Section 5.10 also contained in the demos 
folder of the program. 


A. 1.7 Looping over ranges of parameters 

In order to do the numerical integration for a whole set of numerical points, 
the multinumerics script which was present in version 2 has become obsolete 
in the loop setup. For this purpose, the user only needs to specify numerical 
values for the kinematic invariants in kinem. input, where each line dehnes a 
new kinematic point. 


A.2 General setup 


The structure of the directory general has changed only slightly in version 3 
of the program. The command to launch SecDec is secdec -g similar to the 
loop case. Templates for the input hies can be generated by secdec -prep -g. 

The command secdec -g -p <param. input> -m <math.m> will calculate the 
integral using the parameters specihed in param. input. To evaluate several 
points one can add the additional option -k <kinem. input> to the above 
command. 

The commands 

secdec -g -algebraic and 

secdec -g -numerics will also work in the general setup. 

• param. input: (text hie) 

In this hie the user needs to specify a name for the functions to be evaluated, 
the desired order in e, the parameters for the numerical integration, and he 
can specify further options. The format is similar to param. input in the 
loop integral case, except that by default all parameters occurring in the 
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function are specified in the file param. input and no kinem. input file is 
needed. However, latter can be used to facilitate parameter scans, following 
the syntax of the hie multiparamf ile of SecDec version 2. 

• math. m: (Mathematica syntax) 

Contains the dehnition of the integrand and further options. 
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